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Due to the overbinding that is inherent in existing local ap- 
proximations to the density-functional formalism, certain re- 
action energies have not been accessible. Since the generalized 
gradient approximation significantly decreases the overbind- 
ing, prospects for density-functional-based reaction dynamics 
are promising. Results on the generalized-gradient based de- 
termination of the CH3-CH4 hydrogen exchange reaction are 
presented. Including all Born-Oppenheimer effects an energy 
barrier of 9.5 kcal/Mole is found which is a very significant 
improvement over the local-density approximation. 



I. INTRODUCTION 

Over the past two or three decades quantum- 
mechanical methods based on the Hohenberg-Kohn 
density-functional theory M have steadily evolved in 
scope, complexity, numerical precision and intrinsic ac- 
curacy. A representative albeit incomplete list of some 
of the successes related to gaussian-orbital density- 
functional applications include the accurate prediction 
of molecular geometries pLpl and the accurate determi- 
nation of molecular vibrational frequencies. P-pf Specif- 
ically, in Ref . M , the local density approximation (LDA) 
was used to accurately calculate geometries of H 2 , C2, 
C 2 H, C 2 H 2 , C 2 H 4 , C 2 H 6 , C 6 H 6 , CH 3 C 2 H and CH 4 and 
in Ref. H, Andzelm et al showed that the LDA would 
lead to accurate geometries and vibrational modes for ap- 
proximately ten other well characterized molecules in the 
same size regime. Handy et al [El performed exceedingly 
careful calculations on the benzene molecule and showed 
that density-functional based frequencies within a local 
approximation are particularity encouraging since they 
are more accurate than the MP2 frequencies. More re- 
cently independent calculations by Quong, Pederson and 
Feldman || and Wang, Ho and Wang M have shown 
that the local-density-approximation reproduces the ex- 
perimentally observed vibrational energies of the Cgo 
molecule to approximately 15-30 cm -1 . 

Another interesting success of the density-functional 
theory that is relevent to this paper is the work by 
Salahub et al [Mis 



and more recently Pederson et al |2j] 
and Handy et al [9| on the methylene molecule. While the 
local-density-approximation underestimates the singlet- 



triplet energy splitting [|2||7j|9|], it leads to the correct 
ordering and accurately describes the linear-bent en- 
ergy barrier of 5.3 kcal/mole 0]. Another example of 
an energy barrier that is accurately determined within 
the LDA is the eclipsed-staggered energy barrier of 2.8 
kcal/mole P,fo) in ethane. These barriers are qualita- 
tively similar since the bondlengths and connectivities 
are essentially unperturbed during the distortion. How- 
ever, in a reaction where bonds are broken and reformed, 
the overbinding in LDA may seriously limit the accuracy 
of the reaction surface. For example, using a highly con- 
verged basis set, consisting of a total of 65 even-tempered 
gaussian functions on each carbon and hydrogen atom, 
we have found that atomization energies of nine hydro- 
carbon molecules are overestimated by approximately 0.7 
eV per bond within the local density approximation. pd| ] 
Since a typical reaction energy is approximately 0.5-1.0 
eV, the need for better accuracy is readily apparent. Re- 
cently, Perdew and Wang |1^| and Becke [[13[ have de- 
veloped the generalized gradient approximation (GGA) 
which has lead to vast improvements in describing the 
bonding associated with light atoms (at least). |p|,pdjl4 
For example, using geometries and densities from the 
local-density-approximation, we find the average error 
in GGA bond energies to be 0.1 eV for the hydrocar- 
bon molecules discussed in the first paragraph. fp|pd|] 
These energies have been compared by properly account- 
ing for effects due to zero-point energies. Also, Salahub 
et al H and Handy et al M have shown that the methy- 
lene singlet-triplet energy splitting is accurately deter- 
mined with gradient corrected energy functionals. Due 
to the improvements realized from the generalized gradi- 
ent approximation, the possibility for obtaining reaction 
energies within a density-functional framework appears 
fruitful and some researchers have begun looking at reac- 
tions within the generalized gradient approximation. For 
example, Hammer et al Jig ] have studied the polarization 
and charge transfer that occurs during H 2 dissociation on 
the aluminum (110) surface. 

Since the ability to calculate reaction surfaces quickly, 
precisely and with chemical accuracy are prerequisite to 
a first-principles understanding of reactive processes and 
catalysis, and since the density-functional formulation is 
intrinsically faster than conventional quantum-chemistry 
shemes, this point is investigated here. The density- 
functional based calculations presented here are numcr- 



ically precise, converged with respect to basis set, and 
include relaxation of the core electrons. In this letter re- 
sults on the well-understood hydrogen exchange reaction 
that occurs between a methyl [CH3] radical and methane 
[CH4] molecule are discussed. The experimental bar- 
rier is well characterized and is known to be about 14 
kcal/mole. pq] This reaction barrier has been studied 
within several quantum-chemical techniques by several 
researchers as well. Jl7|-|l9[ 

Many books on the subject of reaction dynamics ex- 
ist p(| and a recent discussion of how to calculate reac- 
tion dynamics from a reactive surface is due to Truhlar 
and Gordon |21| . Here, the interest is in determining 
how accurate a density-functional based reaction surface 
may be and assessing the merits of the generalized gradi- 
ent approximation as compared to the local density ap- 
proximation. With respect to algorthmic issues, even 
in the favorable case where intuition identifies a small 
set (M) of active nuclear coordinates participating in a 
chemical reaction, the quantum-mechanical determina- 
tion of the pathway over a reaction barrier requires ap- 
proximately L M (with L=5-10) times as many electronic- 
structure calculations than are required to determine the 
energy of an N-atom system. Since scaling with system 
size is already a problem for density- functional electronic- 
structure-based geometrical optimizations and ab initio 
quantum-chemical methods, the need for intrinsic speed, 
chemical intuition and novel algorithmic approaches to 
such problems is clear. The relative speed and more fa- 
vorable scaling associated with density-functional based 
techniques is very attractive for these studies. 



II. COMPUTATIONAL AND THEORETICAL 
DETAILS 

To perform the calculations 

discussed here, the all-electron, full-potential gaussian- 
orbital cluster code 22-25(1 has been used. As discussed 
in Ref. J22| , the potential is calculated analytically on 
a variational integration mesh which allows for the cal- 
culation of the electronic structure, total energies and 
Pulay-Corrected Hellmann-Feynman [£4|,26| f° rc es with 
any desired numerical precision. More recently, [g| we 
have incorporated the calculation of vibrational modes 
into the cluster codes. This is accomplished by calculat- 
ing the Hellmann-Feynman-Pulay forces at the equilib- 
rium geometry and neighboring points and using a finite- 
differencing algorithm. Except where explicitly stated, 
moderately large contracted gaussian-orbital basis sets, 
with gaussian decay parameters of Huzinaga |27|, have 
been used. On each carbon atom, I have used ten bare 
gaussians (a=0.1146 to 4232.61) |27J] to construct a ba- 
sis set of 8 s-type (3 of which have an r 2 prefactor), 
4x3=12 p-type, and 3x5=15 d-type functions. On each 
hydrogen atom, I have used six bare gaussians p7|] to 



construct a basis set of 3 s-type (1 of which have an 
r 2 prefactor), 1x3=3 p-type, and 1x5=5 d-type func- 
tions. This amounts to a total of 133 single or contracted 
gaussian orbitals on the seven atom complex. This basis 
set has been developed over many years of local-density- 
based simulations on hydrocarbon molecules, diamond, 
fullerene molecules and tubules and is accurate. [0 As 
described in detail below, we have further tested the pre- 
cision of these results by using a huge basis set of a total 
of 404 even-tempered gaussian functions distributed over 
the complex, and by using this basis we determined that 
the basis set discrepancies are exceedingly small. The ge- 
ometrical optimizations that were required for these cal- 
culations utilized the conjugate-gradient algorithm and 
the simulations proceeded until the forces were less than 
0.01 eV/A. All the results discussed here have been per- 
formed with the spin-polarized density functionals. 



III. RESULTS AND DISCUSSION 

Pictured in Fig. [l] is a schematic diagram on how the 
hydrogen exchange reaction proceeds. As discussed in 
the figure caption, the reaction proceeds from configu- 
ration A to configuration C by allowing the hydrogen 
atoms to relax as the two carbon atoms come close to- 
gether. The intrabond hydrogen overcomes a small bar- 
rier at the transition state (Configuration B) and then 
relaxes downhill to configuration C. 

Pictured in the lower panel of Fig. is the en- 
ergy surface, as calculated within the local-density- 
approximation (LDA). The minimum energy as a func- 
tion of C-C separation is shown in Fig. 0. This corre- 
sponds to the energy of the seven-atom complex as it 
travels along the lowest energy reaction path in Fig. g, As 
argued qualitatively above, the LDA hydrogen-exchange 
reaction energy is significantly suppressed due to the ten- 
dency to overestimate bond strengths. The final result is 
that LDA leads to a classical energy barrier of approxi- 
mately 0.7 kcal/mole which is a factor of 20 smaller than 
experiment. The reactive surface has also been deter- 
mined within the framework of the generalized gradient 
approximation of Perdew et al. pip^ l Using the tech- 
niques discussed in Ref. [jll| the GGA reactive surface 
has been determined and is shown in the upper panel of 
Fig. g. The minimum LDA and GGA energies as a func- 
tion of the C-C separation (i.e. all hydrogenic degrees of 
freedom are relaxed) are compared in Fig. H. Comparison 
of the reactive surfaces and pathways show that the LDA 
and GGA lead to qualitatively different results. Most im- 
portantly, the GGA leads to a classical reaction barrier 
of 8.7 kcal/mole which is reasonably close to the experi- 
mentally known barrier of 14 kcal/mole. p6[ Also evident 
is that the LDA predicts a significantly more compact re- 
active complex. Within LDA, the lowest geometry of the 
CH4-CH3 complex is bound by 2.8 kcal/mole and has a 



C-C separation of approximately 3.27 A. In contrast, 
within GGA, the lowest geometry is only bound by 0.7 
kcal/mole and has a C-C separation of approximately 4.1 
A. It is difficult to determine from experiment, whether 
the LDA or GGA result is more accurate for this complex. 
However, in this regime, the interactions are not expected 
to be too different from the CH4-CH4 Van der Waals 
interactions which are well understood. For methane- 
methane interactions, the carbon atoms are separated by 
4.28 A and are bound by 0.3 kcal/mole. p8| For ener- 
gies of this size, the effects of zero-point motion must be 
taken into account to accurately determine the enthalpy 
(rather than energy) of formation. To determine the ef- 
fects due to zero-point motion, the vibrational frequen- 
cies of isolated CH 3 and CH4 molecules and of the bound 
CH3-CH4 complex have been calculated. For methane, 
two triply degenerate vibrational modes of 1231 and 3082 
cm , a doubly degenerate mode at 1473 cm -1 and a 
nondegenerate mode at 2961 cm -1 are found. These 
agree favorably with the experimental modes of 1306, 
3020, 1526 and 2914 cm" 1 respectively. For the methyl- 
radical the vibrational modes within LDA are found to 
be 537, 3037 (nondegenerate), 1345, and 3158 (doubly 
dengerate) cm -1 . By subtracting the zero-point ener- 
gies of both the bound and isolated systems the enthalpy 
of formation for the CH3-CH4 complex is found to be 
closer to 0.4 kcal/mole. The fact that the zero-point en- 
ergy is larger for the complex than for the isolated modes 
can be qualitatively explained by noting that for the iso- 
lated molecules there are a total of twelve zero-energy 
modes. However, when the two systems complex, only 
six zero-energy modes survive and six additional positive 
energy modes emerge. Assuming that the original non- 
zero modes do not significantly soften, an overall increase 
of the zero-point motion results. 

Based on experience gained in many calculations on 
hydrocarbon molecules, diamond and fullerene-based 
systems, the energetics obtained with the original ba- 
sis set are expected to be very accurate. To further 
address this point, I have used a huge set of gaus- 
sians to ensure the accuracy of these calculations. On 
each carbon atom I have placed eleven even-tempered 
s-gaussians with a=0.1 to 5000, nine even-tempered p- 
type gaussian sets with a=0.1 to 574.349, and four even- 
tempered d-type cartesian gaussian sets with a=0.1 to 
2.56857. On each hydrogen atom I have placed seven 
even-tempered s-type gaussians with a=0.08 to 138.889, 
five even-tempered p-type gaussian sets with a=0.08 to 
11.556, and three even-tempered d-type gaussian sets 
with a=0.08 to 0.9615. Since the cartesian polynomials 
in front of d-type gaussians can also lead to a spherical 
function of the form r 2 exp(— ar 2 ), there are actually 15 
and 10 spherical functions on each carbon and hydrogen 
respectively. Taking account of the degeneracies of the 
p and d sets this leads to 404 gaussians in the basis set. 
The calculation of the energy along the reaction path 



with this basis has been repeated and the changes are 
indeed small with the GGA and LDA barriers increasing 
by 0.6 kcal/mole. It is noted that this does not sug- 
gest that a larger basis necessarily increases the barrier 
and is probably more indicative of the maximum change 
expected. For energy scales this small, one should also 
consider that a small change in transition state geome- 
try could partially offset this small increase. Regardless, 
the results of this paragraph show that the basis set used 
in these calculations is sufficiently accurate to quantita- 
tively determine the improvements that are realized when 
a gradient-corrected functional is used in place of a local 
functional. 

As is well known, the simulation of a reaction re- 
quires the calculation of a classical reaction surface 
and the additional ability to treat the nuclei quantum- 
mechanically. |20| Methods for such simulations include 
path-integral techniques or the techniques discussed in 
Ref. |pH| . Here the interest is in qualitatively determin- 
ing the scale of the nuclear quantum mechanical effects 
and in assessing the intrinsic numerical precision that is 
possible. The vibrational modes have been calculated at 
the density-functional-based transition state which is re- 
alized by separating the two C-C atoms by 2.677 A, plac- 
ing a H atom midway between the two carbon atoms 
and then tying off each of the C atoms by 3 hydrogens 
with C-H bond distances of 1.103 A (configuration B in 
Fig. fy). Our geometry is in good agreement with the 
geometry of Musgrave et al. who find a C-C separation 
of 2.72 A and a C-H separation of 1.089 A. Q It is 
well known that at this transition state one expects six 
zero-energy modes, fourteen positive frequencies and one 
imaginary frequency. If the translational and rotational 
modes are projected out of the hessian matrix prior to 
diagonalization, the imaginary and low-energy modes are 
found to be at 883i, 52, 280, 282, 507, 643 and 646 cm" 1 . 
If the rotational and translational degrees of freedom are 
not explicitly removed prior to diagonalization, the zero 
point motion of the complex changes by only 0.12 kcal/ 
mole but the results from the more accurate projected 
treatment of the hessian matrix are used here. 

As noted above and pictured in Fig. |^ the GGA en- 
ergy of formation required to prepare this state is 8.7- 
9.4 kcal/mole depending on whether the energy is mea- 
sured with respect to the separated molecules (8.7) or 
the reactive complex (9.4). Recently, Fan and Ziegler 
and Deng, Fan and Ziegler [ P9|j30| have performed similar 
calculations on the CH3-CH4 hydrogen exchange reaction 
within local and nonlocal density- functional formalisms. 
In good accord with the results of this paper, they find 
a classical barrier of 11.7 kcal/mole with a gradient cor- 
rected functional and 1.9 kcal/mole for a local functional. 
The 1-2 kcal/mole deviation between their results and 
the results here is significantly smaller than the devia- 
tions between different quantum-chemical methods and 
should therefore be considered as good agreement since 



the numerical methods used were quite different. Ziegler 
and Fan use Slater-type orbitals, the frozen-core approx- 
imation, least-square-fits for the solution of Poisson's 
equation, Vosko-Wilks-Nusair LDA, |nj Perdew's nonlo- 
cal correlation, |B2| Becke's non-local exchange, Jl3| and 
Boerrigter's numerical integration p3j scheme. In con- 
trast, in this work I use Gaussian-type orbitals, relaxed 
cores, analytic solution of Poisson's equation, Perdew- 
Zunger LD A, p4[ ] the Perdew 91 non-local exchange and 
correlation |12| and the numerical variational integra- 
tion mesh of Pederson and Jackson. p3] The calcula- 
tions here and those of Ziegler et al. use spin-polarized 
functionals throughout. Given all the numerical and for- 
mal differences, it appears that there is good agreement 
between the philosophically identical density-functional 
based methods. 

The reaction barrier for the CH3-CH4 hydrogen trans- 
fer reaction has been studied with quantum chemical 
techniques as well and is known to be about 14 kcal/mole 
experimentally. |lq| Several researchers have used vari- 
ous quantum-chemistry algorithms to determine the bar- 
rier and have obtained estimates in the range of 17- 
35 kcal/mole. Jl7]-[l9[ While many more calculations on 
other molecular systems are required before drawing 
definitive conclusions, it appears that the generalized- 
gradient approximation to the density-functional theory 
is yielding results that are as accurate as traditional 
quantum-chemical techniques for this reaction. Further 
from the point of view of scaling with system size, the 
density-functional framework is undoubtedly more ad- 
vantageous. Hybrid methods, which employ fast density- 
functional methods for pathway determination and more 
traditional calculations for the final determination of en- 
ergetics might also be useful. 
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FIG. 1. Schematic diagram of reaction pathway followed 
during the hydrogen exchange reaction. The lowest-energy 
pathway between configuration A and C is realized by allow- 
ing the two carbon atoms to come within 2.677 A from one 
another at the transition state (configuration B). The remain- 
ing six hydrogen atoms relax to accomodate the environment 
mandated by the C-C separation and intrabond H location. 
At the transition state configuration (B), there is one imagin- 
ery root and the classical equilibrium geometry is determined 
by allowing the C-C atoms to separate and the active hydro- 
gen to relax off center. 



FIG. 2. The LDA (a) and GGA (b) reaction surfaces. The 
y-axis corresponds to the separation between the two carbon 
atoms and the x-axis corresponds to the placement of the ac- 
tive hydrogen atom with respect to the bond center. In con- 
trast to LDA, which predicts a vanishingly small energy bar- 
rier, the GGA predicts a classical barrier of 8.7-9.4 kcal/mole 
depending on whether the energy is measured with respect to 
isolated or weakly bound reactant states. Also illustrated by 
the wire-mesh plots is the need to determine the energy, in- 
teratomic forces and vibrational modes on a two-dimensional 
Born-Oppenheimer surface of variable convexity which en- 
hances the computational cost by approximately 10 2 over a 
simple geometrical optimization for the same system. 



FIG. 3. The minimum energy of the CH3-CH4 complex as 
a function of the C-C separation. This potential curve was 
generated by allowing all hydrogenic and electronic degrees 
of freedom to relax with the constrained C-C separation and 
was derived by the lowest energy paths between the separated 
reactant states and the transition state. The geometries are 
schematically presented at the top of the figure. In compar- 
ison to GGA, the LDA predicts a very compact and more 
strongly bound reactive complex. As is also shown in Fig. 2, 
the energy of the transition state is significantly improved by 
the addition of GGA. The horizontal line at the minimum of 
the GGA potential represents the excess zero-point motion 
that the complex appropriates upon congregation. Includ- 
ing zero-point effects, the GGA predicts that the complex is 
bound by approximately 0.4 kcal/mole. The insert shows the 
energy along the entire reaction path. 
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